import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps
from astropy.convolution import convolve, Box1DKernel

E = np.loadtxt('E.txt')
a2 = np.loadtxt('0p5_50K.txt')
aa2 = np.loadtxt('0p5_200K.txt')
aaa2 = np.loadtxt('0p5_600K.txt')


kB = 8.61733e-2 #meV/K
occ1 = 1./(np.exp(E/(kB*50))-1) + 1
occ2 = 1./(np.exp(E/(kB*200))-1) + 1
occ3 = 1./(np.exp(E/(kB*600))-1) + 1

box1 = convolve(a2/occ1, Box1DKernel(11))
box2 = convolve(aa2/occ2, Box1DKernel(11))
box3 = convolve(aaa2/occ3, Box1DKernel(11))

area1 = simps(box1[41:298],E[41:298]) # 41,298,398 for 2, 15, 20 meV
area2 = simps(box2[41:298],E[41:298])
area3 = simps(box3[41:298],E[41:298])

plt.plot(E,box1,color='black',lw=2)  #
plt.plot(E,box2/area2*area1,color='blue',ls=':')
plt.plot(E,box3/area3*area1,color='red',lw=2)

plt.xlim(4,12)
plt.ylim(0,0.15)
plt.xticks([4,6,8,10,12])
plt.yticks([0.00,0.05,0.10,0.15])
plt.xlabel('Energy (meV)',fontsize=20)
plt.ylabel(r"$\rm{\chi}''$ (a.u.)",fontsize=20)
plt.legend(['50K','200K','600K'],fontsize=20,frameon=False)
plt.title(r'$\mathbf{Q}$=0,0,0.5 (r.l.u.)',fontsize=20)
plt.tick_params(which='both',direction='in',labelsize=20)
plt.minorticks_on()
plt.savefig('Chi0p5_ave_boxcar_norm.pdf',format='pdf',bbox_inches='tight',dpi=600)
plt.show()
